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NINE DY, WW. ABSTRACT 


Two problems are presented which are linear on two adjacent inter- 

vals but not on their union. These problems are aSsociated with the 
; ; ; ‘ AX + BF, O<t< ; : 

differential equation X = Joy , DF, #<t<T » where X 18S the matrix 
Cole where F is a 2 x 1 matrix, and where A,B,C, and D are 2 x 2 matrices 
of functions of t. ty is a variable, hence the differential equation is 
non-linear. Problems associated with this differential equation are 
called semi-linear. 

In the first problem, a condition is found on ty and F which must 
be satisfied whenever x(T) is to be a maximum with y(T) fixed. In the 
second problem, conditions on F and tj are found which must be satisfied 
for x(T) to be a maximum for a fixed y(T) and for a fixed x(t ). A numer- 
ical routine is developed which yields successive approximations to the 
maximum value of x(T). 


The basic theory of the methods is presented, and the problems are 


developed in the context of optimum control. 
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i.  Tmereduet Lon. 

Nonlinearities make the solution of optimum control problems more 
difficult. In general, different techniques must be developed for each 
new nonlinear problem. 

The particular problems to be discussed are those in which a set of 
functions x, of a variable t is related to a second set of functions f. 
of t through a differential equation which is linear on each of the in- 
tervals (0,t,) and (t,,T). If t,; 1s a variable, the differential equa-~ 
tion, though linear on each of the intervals (0,t,) and (t,.T), is 
non-linear on the interval (0,T). Such a differential equation will be 
called a semi-linear equation; problems in which a semi-linear differen- 
tial equation occurs will be called semi-linear prébtems. Such problems 
arise in the theory of optimum control; the problems to be discussed 
will be in the terminology of optimum control theory. 


The specific semi-linear problems to be considered arise from the 


AX +, BF, O€t<t, 


CX + DF ta ncoadeuceal: ? where K,A,B,C,D;, and F 
meet 


differential equation X = 


are described as follows: xX is the matrix o whose elements are con- 
tinuous functions of t; X is the matrix i) of derivatives of x and y 
with respect to t. A,B,C, and D are 2 x 2 matrices whose elements are 
functions of t which are piecewise continuous and bounded on the inter- 
val (0,T). In addition it is necessary that B and D be nonsingular. F 
is the matrix (_) of functions of t; associated with F there is some 
constraint, specified in each problem. F is called the matrix of con- 
trol variables, and those matrices F meeting the specified constraints 
are called allowable. 


Associated with the above differential equation are certain boundary 


conditions which X(t) is to satisfy. Generally the initial point is 


given. In addition, one or more components of X may be specified at some 
value of t on the interval (0,T). Curves satisfying the given boundary 
conditions and for which F is allowable, are called admissible. 

In the terms that have been defined above, and for the above differ- 
ential equation, the semi~linear problems to be solved are the following: 
(1) Find conditions on F and on ty which must be satisfied if x(T) is to 
be a maximum for a given y(T), where t, must be determined and where x(t,) 
and y(t,) are unspecified. (2) Add the constraint thst x(t) is to be 
Specified, for a value of t, to be determined. Then (a) find conditions 
on F and ty which must be satisfied for a curve to be admissible; (b) find 
additional conditions on F and on t, such that an admissible curve will 
yield the maximum x(T); and, (c) develop a method of successive approxima- 
tions which will give values of x(T) converging to this maximum. 

The following results are obtained: In the first problem, a neces- 
Sary condition at ty is derived. In the second problem, conditions for 
a maximum are derived, and a numerical method is given for obtaining this 
maximum. These results are another step in the solution of nonlinear prob- 
lems of optimum control. 

I wish to thank Professor Feulkner for the encouragement and help he 
has given me and for the permission to use class notes of his course, 


Differential Equations of Optimum Control. 


II. General terminology and some sufficient conditions for the linear 
problem. 


In this section we define the linear problem and explain the termi- 
nology to be used. We then derive some sufficient conditions for the 
linear problem. 

Let us consider the differential equation 

(1) X= AX + BF 
where X is the n x 1 matrix (x;); where A and B are the n x n matrices 
(a; 3) and (Dia)s respectively, and where F is the n x 1 matrix (f.). We 
will always assume that B is nonsingular, although the results apply to 
many problems where B is singular. The ae the oral and the f; are 
functions of t that are bounded and piecewise continuous on the inter- 
val (0, T). The a;.; 


J 


tions satisfying a given constraint but otherwise unspecified. The f. 


and the Dis are givenifunctions; the f. are func- 


are often called control variables and the Xs state, or dependent 
variables. 

Associated with the differential equation (1) are certain boundary 
conditions which we want X(t) to satisfy. Generally the initial point 
is given. In addition we may specify that certain elements of X take on 
stated values at various fixed points of the interval (0, T). Note that 
the differential equation is linear, since A and B are functions of t 
ayn le 

If we multiply both sides of equation (1) on the left by the 1 x n 
matrix KE (the transpose of the matrix K), whose elements kt are func- 
tions which have not yet been specified, and integrate from t=0 to t=T, 
we get 


T iT a 
(2) \ KIX a = \ Tax dt + \«TBr dias 
9 


a) Ox 


Integrating the left side by parts and rearranging terms, we get 
1h Al T 
(3) Kl x = cs + Kk A)X dt + \ «TBE dt. 
0 0 


Now let us choose ki = . KT and take the transpose; we get 


(4) coy are 


This is the adjoint equation associated with equation (1). Choosing K 


as a solution to the adjoint equation and substituting into equation (3), 


we get 
T et 

(5) KXx| = Ss KI BF dt. 
0 0 


: 1 ole ; 1 . 
We can choose a solution K to the adjoint equation so that K Top) is 
the matrix (1 0 O ... 0). If we substitute this solution into equa- 


tion (5) and rearrange terms, we get a solution for x,(T); namely 


ft 
(6) x, (T) = KT (0)x00) + \ et BF dt. 
0 


In the same way we chose Ki, we can choose ae K, aaa her such that 
Keto) = C6); 6; ooaae 64)» where Fir equals zero, if i#j, or one, if 
i=j. Suppose now that we have the n linearly independent column mat- 
rices Ki chosen above. Such a set of n linearly independent matrices 
is called a fundamental set of solutions for equation (4). If we form 
the n x n matrix K whose i'th column is the matrix on then we may 
write every solution to equation (4) in the form KC, where C is an 

n x 1 matrix of constants. Conversely, every product of the form KC, 
being a linear combination of solutions to equation (4), is also a 
solution. Furthermore K is itself a solution. This is shown as fol- 


lows: Take any product Kc, where C is an n x 1 matrix of constants 


and where K is the matrix defined above. This product is a solution 


to equation (4). yee A = alc. But this equation is valid for 
every C and hence for the C whose transpose is (1 0 O...0). Substitu- 
o 
ting this C into the above equation, we see that the first column of ke 
is the same as the first column of ah But we could have equally well 
chosen the itth element of C as one with the others zero; we hence could 
have found that the i'th columns of both sides were equal, for isl, 2, ..., 
ne —o ~AlX, i.e. AC is a solution to the adjoint equation, which 
is what we wanted to show. Furthermore, if 6 Kis a small variation in 
o o e ° 
K , and ?ift & fis the corresponding variation in Kes then K + § K = 
A iC +6 me , 1.e. 6X = -A'8 7, since ae KX. Hence the varia- 
tions in K are related to the variations in K by the adjoint equation. 
We have shown that solutions to the adjoint equation can be used to 
find solutions to the equation X = AX + BF when F is known. The next 
problem we are concerned with is that of finding F so that for a given 
X(O) and a given T, a specified component of X(T) is a maximum. F can 
be regarded as a column matrix of forcing functions; these forcing func- 
tions are the components of the forcing function vector. There may be 
one of several types of constraints on F. [In rocket thrust scheduling, 
for example, the acceleration possible at any one time is limited by a 
function of the mass of fuel on board at that time. If we call this 
function g(t), the corresponding constraint on F is that lrle» G).., 
where |F| is the Square root of the sum of the squares of the elements 
of F. Another type of constraint on F is [ealé [ail where a, is a 
given function of t. Problems of this second kind are called bang- 
bang control problems. {i} i In every problem it 1s necessary to state 
the constraint on F; forcing functions which satisfy the stated con- 
straints will be called allowable. Solution curves to equation (1) 


for which F is allowable will be called allowable curves; allowable 


curves satisfying the differential equation and satisfying the given 
boundary conditions will be called admissible. 

We want a method for choosing an allowable F such that X(0) is a 
given point and such that a specified element of X(T) is a maximum. A 
most important principle enables us to state conditions which F must 
satisfy whenever the desired maximization takes place, namely Pontry- 
agin's Maximum Principle: The control variables must be chosen from the 
set of allowable controls so as to maximize a scalar product of some 
solution to the adjoint equation and the forcing function vector at every 
time tu[ZPR we will use this principle extensively in the following 
pages. 

Consider the following problem and see how the maximum principle 
may apply to it. Suppose that we want a curve beginning at some fixed 
point at time t = 0 on which some element of X at t=T, say x,(T), is a 
maximum and such that on the curve these three conditions hold: First, 
F is allowable. seca the curve satisfies the differential equation 
X = AX + BF for all values of t between t=0 and t=T where T> 0 is given. 


Thee , ml) = Xin for 1©2, ... ,; Ny WHeBe x. 


it are given, and where 


1<Nsn, and where X47 6T), eG x, (T) are unspecified. N=l means that 
no values of the x; are given at t=T; N=n means that all values except 
X, are given at t=T. A proof given by Faulkner (33 proves that the fol- 
lowing hypotheses are sufficient for a curve C* to yield a maximum x, (T). 
Suppose that we have found a curve C* with forcing function F* and have 
found at the same time a solution K* to the adjoint equation which to- 
gether satisfy the following hypotheses: 

Hl. C* is admissible: it begins at X, and ends on the manifold de- 


fined by X(T) = Xoms co% 9 Xp CT) = x and F* is allowable. 


NT’ 


10 


H2. F* maximizes K*=BF for all allowable F, i.e. K*=BF*> K’ {BF for 
all allowable F. 
H3. k}(T)>0, and k; (T) = 0 for i>N. No restriction is put on 
* a 
k; (T) for L = 25 eee yN. 
THEOREM: C” furnishes the desired maximum of x,(T). 


re 1 
Proof: We have shown that KI x = \x'BF dt for every solution K to the 
0 0 


adjoint equation. Hence for the particular solution K* and the matrix 
X* we have, on the path c™, 
* * x * a eT x 
C7) k1(T)x, (CT) + Ky (1) X on PP Ky OD) Xa = (K X dm 
wail 


ii 
* 
= (K x”) + \x ThF dt. 
ss 0 


Consider any other admissible path C' with F = F'. For this path and 


for K = K* we have 


: * xT x kp 4 
(8) k](T)x1(T) +... + ky(TXyp = (KX Dg + (K BF at. 
0 


Subtracting equation (8) from equation (7), we get 


T 

t . * * 

(9) k} (x, -x1) 7 * \¢K TRF -K BF Jato 0. 
0 


Hence, since ky CE) > 0; x1 (T)> ms at t = T. Hence the given hypotheses 


are indeed sufficient to give a maximum x, (T). Note that since kK’; (T) = 
0, for j=N+1, ... , n, the unspecified elements of Xx" (T) play no part in 
the solution. 

Having developed a useful sufficiency condition, let us now consid- 
er the case where X is the matrix ()» where A and B are 2 xX 2 matrices, 
where F is the matrix (@). and where we want to maximize x(T) subject to 


the following constraints: First, (£1)? ~ (ie = 1. Second, TO is 


Ll 


fixed. Third, solution curves must start at _ and end on the line y(T) = 

Ye» where yr is given. For this problem the admissible curves are allow- 
sd 

able curves satisfying the third constraint. 


Let us choose the solutions K} and Ke to the adjoint equation having 


K1(T) =(}) and K2(T) = GP Using these solutions in equation (6), we 


get 
(10) x(T) = (KD), + \K "BF ae 
0 
and 
2 a 2 
(11) y(T) = (K iy + \K BF dt. 
0 


Since the given constraint on F is that eae, * = 1, we may choose 
fy = cos 8, f. = sin ®, where @4is a function of t. With this substitu- 
tion, a variation in 98 will cause a variation in F which will in turn 
cause a variation in x and in y. Let us call the variation in 8, 6®@ and 
the variations in x and in y at t=T, 6x and Sy, respectively. We then 


get, from equations (10) and (11), 


T 
iM 
(12) bx = \K Be 6@ at 
and 
(13) bby = Jx7 "BE 60 at, 
0 


for sufficiently small values of 66. Let us now choose a special varia- 


oF. 


: = Las opt ye 
tion of the form 6 = ejK Bsn Ct e,k BO where e, and e, are con- 


Stants yet to be chosen. Substituting this variation in equation (12), 


We Pet T 
(14) px = Ck pee ft OF 
= ae) es BSo + eoK Bg) dt 
0 


12 


and a similar expression for §y. These can be simplified if we let 


T 
tJ = \ TRE) (Kp Fat. 
0 


Using this substitution we get 


11 12 
use 6x = e,! 4: en! 
and 
a 24 22 
(16) Cy topic l 
; 12 ag : Vy 
Note first that I = I. Note next the consequences if I = 0. If 
2? = 0, then K2Tpe® F must be identically zero; hence t2! and By are 





9 6 


both identically zero; hence y(T) is not affected by a change in the con- 
trol variable 8. We describe such a situation by saying that the curve 


furnishes a stationary value for y(T). We will not want y to have a 


22 


Stationary value, and we shall henceforth assume for our curve that I ‘# 


Tp oF 


67 18 not identically zero on the inter- 


0, or, equivalently, that K2 
val (0,3). 

Since T is fixed, equations (15) and (16) give the total differen- 
tials of x and of y. If we replace&y by y¢z - y(T), we can generally 


find values for ey and e, which give a variation which will make the re- 


2 


sulting curve admissible. 


Suppose now that we have an admissible curve and that we want to 


, wed ia] : : ee: 
find conditions on the I y such that x(T) is a maximum. For an admissi- 


ble curve we have 


wd 
(17) 6a = ey 0 + a1 


and 


ear +e t2*. 


Gls) 0 D 


Equations (17) and (18) can be solved for §x»0 only if the rank of 


irs 


Aileen? iil <i 
if L iE E 
ia is equal to the rank of 4 - But the rank of the first 
6x 0 0 
matrix is the same as the rank of - The rank of this last 
0 12! 122 


Zi 172 


matrix is one greater than the rank of the matrix (I ); hence equa- 


tions (17) and (18) cannot be solved for 4x whenever the rank of 


11 ti2 1 22 


is equal to the rank of (I I). But this will Dee@mwevonly 


I 


5] is zero. Hence if the admissible curve yields 


ri5| = O. 


when the determinant 








a maximum for x(T), then 
Let us consider some of the consequences of this condition. If the 


determinant 





ri5| = 0, then the matrix (I74) has a zero eigenvalue. Let 
= 


us choose constants Cy and Cc, such that (: 
2 


is an eigenvector of (177) 


corresponding to this zero eigenvalue. Then let us consider 


(19) scope Tage, + cK pOFy2 at. 
08 2 38 
Clearly J20. But J = Como + 2cyep1 1 Coa” 
21 122 4 
=(cies) : 
2 12! 122 e, 
my 1) 
But 1S an elgenvector corresponding to the zero eigenvalue of (I ~), 
Cc 
2 


i 2T : : 
hence J = 0. Therefore c_K Ea Cak pal =o. This equation 
1 o 8 2 
appears frequently in the literature on calculus of variations - it is 
known as the Euler equation. Curves whereon the Euler equation is Ssat- 

aunt: 1 x 1 ja 
isfiled are frequently called extremals. Hence K = c}K + coK’ 1s a 


this is the definition given by Bolza (iy) . Another definition is 
given by Bliss [5} 


14. 


solution to the adjoint equation which gives an extremal. 


5 


* 
We saw that if c,>0 and if @ = @* maximizes K ~BF then x). sea 


Maximum. The condition KT as = 0 is a necesSary condition for x(T) 
to be a ttaximum. We can and will choose C1? O. We will also see that 
it is necessary that ® maximize K BF; this is the Weierstrass condition. 
In this problem the Weierstrass condition is necessary and sufficient 
for an admissible curve to furnish the desired maximun. 

Let us assume that the rank of (14) is at least one for all curves 
which arise: this is called normality for the problem. Suppose now that 
we have an admissible, normal curve which furnishes the desired maximum. 


On it the Euler equation must be satisfied, and the matrix qa has a 
zero eigenvalue with a corresponding eigenvector (22) s we will choose 
c;>0. Then the F determined by the Euler equation maximizes the pro- 
duct K lpr, as a function of 8, over the entire interval from t = 0 to 
t=T. The proof is as follows: 

Suppose that 8, (t) is the argument of F for some admissible curve 
satisfying the Euler equation but that 9,(t) does not maximize the pro- 
duct KT pr (@) on some subinterval (tz sto) of (0,T). Then there is some 
other 6, say 65(t), such that KT BF(®,)> KBF(8,) on the interval (tj, 


to). | Let O= O65 on the interval (ty,t} radii), for ty + dt, <t,; and 


let 66 = eqk BSH + ek? Boe elsewhere on the interval (0,T) . 


Then, since both paths are admissible, 


igi 11 12 
(20) dx = K B(F(@,) - F(®1)J dt, + Ive, + Ie, 


and 


(21) 0 


2T 21 22 
K° BE F(@5) - F(8,)) dt, +I e, + I°’e,. 


If we now multiply equation (20) by c, and Gal) *by Co and add, remember- 


ing that KT = i + 2 ale we get 


15 


(22) c,dx = KB |rce,) - £(81)} at, 


since(.)) is an eigenvector of (144) corresponding to its zero elgen- 
value. Since the right-hand side of equation (22) and the constant Cc) 
are both positive, the dx of equation (22) is positive. Furthermore it 
is possible to satisfy equations (20) and (21) by a set of e; which will 
give a positive dx and at the same time keep dy = 0. To see this, mul- 
tiply equation (20) by Cy and call the new equation (20'); multiply 
equation (21) by cy and call the new equation (21). Add equation @o ) 


to equation (21') and call the resulting equation (23). Looking now at 


: : mT y= 2 
equations (23) and (20), remembering that K ~ = c,K | #80 OK » we have 
* 
(23) e,dx = KB [F(05) = F(07)) at, 
and 
(20) Ov= KB [ F085) - F(8,)] dt, + Rapes, + i ae 


But, by assumption, KB [F (82) ~ F(81)] >0, and 22 x 0, hence we can 
solve these two equations for aT and en: Hence it is possible to satis- 
fy equations (20) and (21) by a set of e; which will give a variation 
which will in turn yield a larger value of x(T). Hence if @(t) is such 
that the product K* TBF(e) 1S not a maximum on every subinterval of (0,T), 
then x(T) will not be a maximum. It is worthy of note that the form of 
the conditions for an extremum is independent of the particular problen. 

The Weierstrass-Erdmann corner condition is an immediate consequence 
of the above condition. For suppose that we have a curve whereon the 
hypotheses of the above theorem are Satisfied. Suppose further that tj 
is a point at which the control variable @(t) is discontinuous. Then 


t 


Cec KIT 2T 3 = 0 f: : ; . 
1 C5K* ") BF (8) a = 0. Proof: This must be so; otherwise the 


condition just found would not be satisfied in some neighborhood of ti> 


16 


for either t> ty or t<t,- The condition that K BF (@) “ = K* I pp 
1 


(6) pps is commonly known as the Welerstrass-Erdmann Corner Condition. 
In this section we have introduced several important concepts as 

they pertain to a linear problem discussed in the following gections. 
In the next section we consider a combination of linear systems, 


the combination being non-linear. 


ty 


17 


III. The Weierstrass~Erdmann corner condition for a more general problem. 
In this section we consider a problem having a corner at a variable 
time. We develop the Weierstrass-Erdmann corner condition for this par- 
ticular problem. 
Let us consider the differential equation 
(24) X= A'X +BF 
where A' is A, for 0<t< ty> and C, for ti< t<—T, where B is B; £0r 
O<t<t,, and D, for ty<t<T xX is the natrix(") ; F is the matrix 
(E2); A,B,C, and D are 2 x 2 matrices of functions of t which are piece- 
wise continuous and bounded on their respective intervals; B and D are 
non-singular. T is fixed; t, is a variable to be determined. The con- 
Straint on F is as before, namely (pr + (e)° = 1. In addition, x 
and y must be continuous at t,;. For this problem the admissible curves 
are the allowable curves starting at X, and ending on the line y(T) = Yer 
Note that this problem is non-linear on the interval (0,T), since 
A’ and B’ are functions of t, as well as of t. The problem is, however, 
linear on each of the two sub-intervals (0 ,t4) and (t,,T)- Hence we call 
this a semi-~linear problen. 
We can rewrite equation (24) as 
: AX + BF, OSt<ty 
(25) X = 
CX + DF, t7<t<T. 
The adjoint equation is 


-ATK, O< ticeer, 


(26) K - 
=0°K, t, ema 


. il 

We choose, aS particular solutions to K = -clk, the solutions K1(T) = (") 
0) ° 

and K2(T) -(°). For ti< t<T, then, each K? is a function of T and of t. 


Since T is to remain fixed, however, we shall suppress it wherever it 


18 


occurs in the Ki and shall consider the Ki as functions of t only, on the 
interval (t,,T)- Let us define a set of particular solutions to the 
equation K = alxK by Kae, = aCe for 1 ="1,2. For t°on’the inter~ 
val (0,t,), then, the Ki are functions of the variables ty and t where, 
aS above, we Suppress #. 

Before continuing further let us make one substitution. When con- 
venient, we shall use E to mean B, in the interval (0,t,) and D, in the 
interval (t,,T). (Note that E and Bare the same matrix). Using this 


convention, we get 


a 
(27) re) K "(ty ,0)X(0) + ("be dt 

0 

and 2 
(28) y(T) = K "(ty ,0)X(0) + \ x7 Tae dt. 

0 


2 
Since the constraint on F is that ne + (£5) = 1, we can replace fy 
by cos © and f5 by sin 8, where 9 is a function of t. Since we are con-~ 


sidering both @ and tj as variable, we get, from equations (27) and (28), 


._ &x(T) 
(29) dx(T) = 2 tn + 6x(T) 
and 
_ dy(T) 
(30) dy (T) a dt, + 6y(T) 


where §x(T) and $y(T) are the variations in x(T) and in y(T) due to 
variations in the control variable 98. 
§x(T) and éy(T) are given by the following equations for sufficient-~ 


ly small values of 98: 


4p 
(31) &x(T) = (ct 3 68 dt 
0 


and 


19 


T 
(32) 6 y(T) = Nis 68 dt 


Let us choose the special variation 60 = e,k ES + e)k ES where 
ey} and @5 are constants which remain to be specified. Then equation (31) 
becomes 


ab 
Fag 
(33) ox(T) = VK BRE) (ek ERE + ek ESE) at 
0 


with a similar expression for 6y(T). Adopting the substitution 


T 
cays s Vote RE) (ESE at, 
0 
we get 
(35) X(T) Fe,1"? + be 
and 


21 212 
e,l + eol F 


(36) 6 y(T) 
We will use equations (35) and (36) later. 

The variations in x(T) and in y(T) due to a variation in t, are 
more complicated. Rewriting equations (27) and (28) to show the way in 


which t, enters, we get 








e) aT 
C37) x(T) = “T(t, 50) X(0) + Jere, BEE + auesy dt 
| 
and ty T 
(38) y(T) = 2F (t, ,0)X(0) + VK Te, ey BRS + \x2"ceype ae. 
0 ti 
Hence a 
ie 
Ox) at, = 29K _¢t,,0x(oy at, + ( 94 
39 — = 3 
(39) ore a 1 St (t,,t)BFdt jdt 
0 


1T 1T 
+| (K BF) -(K DF) dt 
L ti ee i 


and a similar expression for oo » the only difference being that 
1 
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2 ; ' : 
each K!T is replaced by K To avoid problems in notation let us denote 
LT ae : 1T , 
by 6K” the variation in K'~ caused by varying ty by a small amount dt.; 
ee: it _ exit 
i.e. to a first-order approximation 6K ae one With this nota- 
al 
tion, equation (39) becomes 
ys! 
quo) OX at, = ox? t, 0)x(0) (oxi ,t)BF dt 
ot, t s il 


0 
1 1 
+ ((K BE), -(K "FD, J dt, 


Since Ki is not a function of t, on the interval 51), 6K = 0 on that 
interval. For t on the interval (O,t,); aK: Satisfies the adjoint equa- 
tion, and hence ik: = -A'6 xl on (0,ty)- Keeping this in mind, let us 

consider the integral in equation (40). From equation (25), BF = X - AX. 


With this substitution, the integral in equation (40) becomes 
fl 1 . 1T 
(41) (ox * (ede ae = \ox (t,,t)AX dt. 


Integrating the first integral by parts, we get 


ty, ty 


1T i i 
(42) 6K (t,t X(t) (ok 48x dex at. 





0 0 
But ox! satisfies the adjoint equation, hence the integral in (42) is 
zero. Hence (42) reduces to BK (t,t) X(t 1) - BK (t, ,0)X(0). But the 
variation in au due to a variation in ty» evaluated at tj; is equal to 
Lk Te, DA - eT, al dt, > by equation (25). Furthermore, X(t) is con- 
tinuous, so that X(t, _) = X(t ,,)- Hence (42) is equal to 
[kT a - eT, ral X(t 4) dt, - 6K'*(t, ,0)x(0), and equation (40) 


reduces to 


2x(T) eee 
SS a a [ax(t,) + Br(t,_)) dty 


21 


pi 
~ KE gp (CXC) + UE Ly de. 


But, by equation (25), this says that 


IF bd ith e 
(44) <a = |x (t,_)X(t,_) - K (tice, J dt, - 
Similarly, 
(45) — dt; = [« (ty X(t, _) - K Ct ket, dt.- 


Hence equations (29) and (30) give, as the total variations in x(T) and 





y(T); 
© 
¥ Ti oF (xT; 1+ 
(46) dx (T) - | K E-<A— 58 dt - (K x) | dt, 
0 = 
d 
an 7 - 
" 2T. aF 2nite i 
0 ty 
: ie ai 2T, sak ‘ 
Choose the special variation 56 = (e,;K + e, EO - Then, using 
equations (35) and (36), we get 
Le 
: Lt 
(48) d(T) = 1) eet” SS Ce dt, 
ti 
21 22 oT., |r 1+ 
(wy ay(T) = 174e, + 172e, - cM) | at, 
c 
; i 


Observe that if t,. is fixed we have the problem discussed earlier. 


1 
Let us assume, then, that admissible curves exist and that our prob- 
lem has a solution. If admissible curves exist, it is possible to find 


€,, e, and dt, from equations (48) and (49) and hence to get a variation 


2 
$8 such that the curve obtained by replacing 8 by @ + §@ is admissible. 
Suppose that we have performed these calculations and have obtained an 


admissible curve. For this curve, dy(T) in equation (49) is zero. Be- 


fore continuing further, let us assume that Tee 7 0. This assumption 


2. 


is enough to insure normality in this problem. Having, then, an admis-~ 
sible normal curve, we can use the results obtained after equation (19). 
We showed that on an admissible normal curve for which x (T) is a 
maximum, the matrix (tJ) has one zero eigenvalue and that the components 
c, and cy of the eigenvector can be chosen so that c, is positive. We 
also showed that one solution to the adjoint equation which gives an ad- 


- + Sh eee Hence let us choose Cy and C5 as 


missible extremal is K” = cK 
components of an eigenvector of the matrix (r1J) corresponding to the 


zero eigenvalue, with c,;>0O. Multiplying equations (48) and (49) by 
1 





these c; and adding, remembering that the Kl are continuous at th: gives 
Cis 
kT ‘ 
C30) c,dx = K (t,) ee ce : 
1 
This leads to the following theorem: * 
1+ 
THEOREM: For x(T) to be a maximum, the quantity Keay) [kct) ] 
fe 
le 


must equal zero. 

Proof: If the above quantity is positive, any dt,? 0 will give a posi- 
tive dx(T) and hence a larger x(T). If the above quantity is negative, 
any dt,;<0 will give a positive dx(T). This is the Weierstrass-Erdmann 
corner condition for this problem. 

The next section uses the problem we have been studying for back- 
ground and considers the essentially different problem obtained by intro-~ 
ducing the constraint that x(t ,) has some fixed value. We develop a 
numerical routine for determining the solution by a method of successive 


approximations. 
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Iv. A numerical routine for determining the maximum x(T) for a fixed 
value of x(t) 


In this section we consider the problem states as before but with 
the different constraint that x(t,) has some fixed value. We make up a 
curve consisting of two arcs, each of which is an extremal. We use the 
method of variation of extremals on these arcs to drive the resulting 
curve to admissibility and in a gradient technique to determine the curve 
on which x(T) is a maximum. 

Let us again consider equation (24), namely 

(24) X = A'X + B'F 
with the same conditions as in the beginning of Section III but with a 
different constraint on ty> namely that x(t,) = X,, where x} is given. 

This problem, like the one in the preceding section, is semilinear. 
It is not linear on the entire interval (0,T), since A’ and B' are 
functions of t, as well as of t; it is, however, linear on each of the 
two sub-intervals (0,t,) and (t,,T)- 

The differential equation we have been working with is 

AX + BF, O0<t<t 


1 
<t< F 


G25.) X 


CX + DF, ty 


For this problem the admissible curves are allowable curves satisfying 
the above differential equation and such that X(@ = Xo» x(t,) = X)> 


and y(T) = Yeo where x_ and Ye are given constants. The adjoint equation 


1 
for equation (25) is 
~Alk, O<t<t, 


(26) K = Fc 
-C’K, t,< ts T. 


1 


We choose solutions K~ and K* to the adjoint equation such that Kl(r) = 


1 0 - - 
(>) and KT) = (’) and such that &*:ate continuous at t = ti- We further 


define the 2 x 2 matrix K as before; its first column is K (+) and its 
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second column is K2(t). Note that W is a function of the variable t for 
t,;<t<T and of both the variables t, and t for O<t<t}- 

Multiplying equation (25) on the left by ae integrating from t = 0 
to t = t,, and using the fact that YW is a solution to the adjoint equa- 


tion leads to 
ty t 
(51) Ce txt) | = VE (ty toBF at. 
0 0 


We have also, from the last section, 


t T 
1 
yy wo axCie= K Tceemoncore ( kK’ (t,,t)BFdt + \x Tce) peat 
0 ty 
Ly ak 
(38) y(T) = K7(£,,0)x(0) + \ KC, ,t)BFat + {7 ceyprat. 
0 t 
1 


The differential equation is linear on each of the intervals. 

We choose F by using the method of variation of extremals on (0 ,t4) 
and each of them (t,,T)- On each of the two arcs, we have the following 
by the theorem of section II. If a. is an admissible arc, if K" is the 


: — : ; * 
solution to the adjoint equation defined by K = ec Kl 1 


2 
1 + CoK » where K 
and K* are the solutions to the adjoint equation defined above, if C1) 0; 
‘ * * — xT * P 
and if, on C , F maximizes the scalar K “EF, then C furnishes a max- 
imum x at the end of the arc, relative to the point at which the arc 


began. If we, then, choose F* to maximize the product Kt 


BF, for 0€t 
<tj, and Kl pp, for t; <t<T, the resulting curve will be made up of an 
extremal from t = O to t = t, and an extremal from t = t, tot = sie 

Once we have an admissible curve made up of two extremal arcs, we 
want some way to vary these curves so as to get a maximum x(T). We use 


a gradient technique, explained later, to choose a set of variations in 


the curve parameters to get a larger x(T). The new curve may not be ad- 
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missible because of second-order effects. Hence, we again drive the 
curve to admissibility, say, by the method of variation of extremals 
mentioned Pe We continue this process until the maximum x(T) is 
obtained. 


The calculations proceed as follows: Since Cy ls positive on each 


1 


* Z ‘ 
arc, we can choose Cc, = 1. Hence K = K” + cK’, where c 18 a constant. 


Unfortunately when a condition on X is given at time tj, K* may not be 


continuous at t)- Hence let us suppose that c 1S not the same in both 


arcs, and let us adopt the following notation: K* = K} + ak, for t on 


1 


(O,t)); Ky = K* = bKZ, for & on (t,,T)> where a, b, and t, are unknowns 


which must be determined so that the curve made up of the arcs is admis- 


*T RF 


sible and yields a maximum x(T). Now consider F". F* maximized K 

on the interval (0,t1)3 on that interval K” is a function of t,>t, and 
oe : 

a. On the interval (t,,T), K is a function of b and t only, hence F* 


is a function of b and t on the second interval. 


With F replaced by F’, equations (51), (37), and (38) become 


Gc 
(52) A(t, t)X(T) ' 2 f A (t,,t)BF*dt 
ty T 
(53) (TY SK" (tgs 0K COD \xlT ce, 2 BF*at q (Te prtae 
0 £1 
ta T 
(54) y(T) = K~'(t,,0)X(0) + WT Ce, tBR*at ‘ aera 
0 ty 


From these equations we want to devise a routine that will, first, give 


Itt is possible that we might obtain an x(T) from a set of curves 
that are admissible under one admissibility criterion but not under 
another. After trying various admissibility criteria, the author de- 
cided upon that of calling a curve admissible if [x(t,)-x,| + ly(T)-y-| 
$10 °. 
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us a set of admissible curves and, second, find an admissible curve where- 
on x(T) has its maximum value. We wart also to find a condition which 
will indicate that no admissible curves exist, if such is indeed the case. 
In this problem, a curve can be characterized by a set of values for ty> 
a, and b. Hence we want a routine to find values of a, b, and ty which 
will determine an admissible curve whereon x(T) is a maximum. 

Let us first get expressions for the total differentials of x(t,)>, 


of x(T), and of y(T). First, of &(t,). Equation (52) can be written 
ty 
T iE * 
(55) A (t,t, )xX(t,) = 7 Ct, .0)X(O) ADEE at. 


Taking differentials of both sides gives 


(56) dAMP (ey .t,)X(t,) + ATE, t pax(t,) = aK (t,50)xC0) + 


t 
1 
T # dé * 
yd (ty,t,)BF (t,_)dt, + ( dx (t,t) BF dt + 
0 


ca ‘Af * 
‘ K(t,,t)Bor dt. 
0 


To get i Crete note that the first t) denotes the end of the inter- 
val (0 ,t4) and that the second is the value that the running variable t 
assumed at the end of the interval. The differential due to the change 
in the first ty is, since AW satisfies the adjoint equation, (8 Gee 
(A-G)dt,; the differential due to the change in the second t, is =A (ty, 
i) Audt Hence axrct eB ) ee Act t.)Cdt.,. and aac t.)mx(t.) = 

1 i” jr Le il 1 aL 1 

- H(t, ,t,)CX(t,)dty. 


1 
a 
Consider next a A(t, ,t)BF dt. By the argument after equation 
0 


(42), this integral is equal to A" (ty st) (AC) X(t, at, ~ dA T(t, ,0)x(0). 


Substituting the above results into equation (56) yields 


va 


(57) (tpt aXe) = WCE, t,dAK(t, dade, + 
a 


Ac" (ty ,t,)BF(ty_)dt, + \ AC (t, st BOF dt 
0 


t, 
=x! (ty ot )X(ty_ if x (t, ,t)BOF dt. 


Calculating the integral in the above equation is somewhat more diffi- 
Ke al oe al 2T ; 
cult. F is chosen to maximize the produce (K + aK )BF. Expanding 


the product Ce + ak?!)p yields the matrix 
1 Min. li, 12 oT iad ) 
( (ic sak’ 7] by, + [ko+ak",, fc teak tI b, 5+ (x? t+ak?/ v,,,, 


Calling the first element hy and the second h, and remembering both that 


2 


F* must maximize the product Ch, ho) F and that F can be written as 


sin 8 


2 9935 x  f=sin 6 : = 
[(hz) + (ho) | - Hence 6F = ~ 66, where 6 tan Ch, /h,)- 


1 
(*". ’); we see that cos @ = hj/ [ cn) 2 + (hy)? ] 2 and that sin 6 = ho/ 


Substituting in the expressions for sin 8 and for cos 8 and finding 66 


in terms of 6h, and 6h5 ylelds 


oot (*" hybh, ~ hy oh, 
hy [ch,) +(ho) 7 


oh, and 6 ho are each the result of variations both in t, and ina. Con- 


sider first the part of the variation Ch, oh ~-h, oh » due to a variation 


2002 
in a. 
oho oh 
a2 Me re le 21,_,22 12 22 7 
1 2 2a 22 12 22 
- [cx +ak )bL5 + (k +ak Yoo] [k by4,tk boy 
= 12 22 12 22. 22 
= a [(k*b, +k**b 54 (kb, +k.) - (x *b, atk °b, 5) (k* oy, 


22 


+k 12tk ia >) - Ge 15 


22 bo] + (Geb, +k2 4b.) kb " 
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21 12 22 
11.12 21, 22 
a Eee) ge Ween 9? 22P21) 


11. 22 21.12 
+k k (by, bo9-by obo) + KK (bp 1 by 9-bo9b1 4) 


eel 1 gad 2.2 
= (kook **-k ok) (by by 9-by 5b5y) 


KI > |). 


where |X| denotes the determinant of K and similarly for [3] . Hence 


the variation due to a variation ina is|K| | B| da. 


The part of the variation (hy Shy-hoohy) due to a variation in t, is 
(59) ae = ve dt, = fcc tak) ?yb, 1 +(e? ak??) bd, 
{cen + abk *)b, 5+ (6k? b+abk*)by 5] - 
fact + ak’“)b,, + (k? +k"), 0] [cox 


tak )b,, + (6° + adk**)b,, ] 


2 atl 
pity ak l? ape an 


“3 
21 
kas ak ~~. MAK tnask22 


For convenience, we shall refer to this variation as VAR(1) henceforth. 


Note that each of the 6k / can be calculated, since 6 A(t, »t,) = (atc!) 


1G eae, and since 6 K= -alg Kon (0,t,)- 


It is possible to simplify, within the integral, the terms KC (te, ,t)B 
(72): as follows: 
1 
1 ee 72 11 12 21 De 
-h, k k bi 19 =. (k .+ak )b, +k +ak boo 


T 
AC (1,08 . 
hi} \Kl? «777 \b,. b (kt beak’ *)b, 1 +(e? b+ak??)b, 5 


21) ae 
-a 
1 





| lk 
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Using the above substitutions, equation (57) simplifies to 


(60) mit, .ti)dx(t,) = AC (ty st Kt) + 
| 
-a | | 
\ ( se rat | da + VAR(1) dt age 
tl [cn,)2+(ny)2] 3/2 


*! is a non-Singular square matrix, hence (Gt) =" exists, and 


° gh -1 
(61) X(t) = K(ty_) dt, + [A et, * 


Sil 
= | IBlda + VAR(1) dt 
B Serr tt t, 
") { iin, )* of : 


4 + (ho) 


dx(t P 
where dx(t,) = = Ay - For future reference, let uS rewrite dx(t,) 


dy (tj) 
as 
(62) dx(t,) = a, dt, + a,5da; 
do x(t1) @ x(t) 
where a = —————— and a = —=———" .__ Thus we have derived an ex- 
le oo poe ar a 


pression for dx(t,) in terms of variables which we can calculate. 
We want next to derive an explicit expression for dx(T). We can 


rewrite equation (53) as 


* 
ae ae 
1 


ee Tol 


(63) x (T) = K' "(ty ,t4)X(t,) + 


whence 


Gey PP asen = dk h(t, separtt.,) i K (ty st aKxct,) 


ily 
Tt * il 
- K (t)DF (t, dt, + ( K "(+ }USF “dt . 
t 
1 


By the discussion preceding equation (57), dk" (t,t) = ak Gagt,) Gata), 


IT, 


1 
and dK y2 tz) X(ty) = -K “Go NORE ae. If we write the differen- 
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ae 


tial for dy (t,) corresponding to equation (62) as dy(t,) = s% | 1 


Aeoda, the product KT (t,t, ax (t)) becomes 


12 


1 12 11 
(k ce +k “a,,)dty + (ko O), + k'“Gs5)da. 


sky 
* 
By an argument similar to that preceding equation (60), noting that F 


is dependent only on b in the interval (tj,T), we get 


(65) \ KT(typor*at = \ (-b) _|B|* [K [Pat yap 
25) uy kn) 2+thy) 2°”? 


Combining the above results yields 


1T * 11 12 
= a (@ & 
(66) dx(T) K (tot) [Cx (t |) +DF (t, Jat, (kG, tk “A, dt, 


T 
+k? ha) +k a) das \-b) [pl * fel@at ab. 
(ch) +Gh,) 57 2 


=) 
Combining terms gives us an equation corresponding to equation (62), 


namely 


dy = 22D, and yy = 2D. 


dy(T) can be derived much as we derived dx(T); the equation corresponding 


to equation (66) is 


a 2T * ‘ 
(68) dy(T) = -K°°(ty,t,) [CX(t,) + DF (t, Jat, 
21 22 21 22s 
+ (k @4,*k Agy dt, + (k @1o+k a 5o)da 
T 
“\ Ini? [X12 at ap, 
2 2973/2 
[cn )7+#¢hy) J 
25) 
which we abbreviate as 
=o ef a : 
(69) dy (T) aptiliiyp te Capi ol 
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We thus have three equations to use to determine corrections in ty» a, 


and b which will lead to a maximum x(T), namely 


(62) dx(t,) = a,,dt, + a, oda 


a5 dt, + ada + a, ,db 


(67) dx(T) 


a..dt. + a..da+a._.db 


(69) oe) 31°°1 * %39 33 


The numerical procedure used for solving for a maximum x(T) is the 
following: Choose values for t)> a, and b. With these values, calcu- 
late x(t ,), x(T), and y(T) and eet G2 starting from X(0) = X,-  Prob- 
ably the curve so determined will not be admissible, i.e. x(t) wlll not 
be xy and y(T) will not be Yeo To obtain a curve which will be admissi- 
ble, solve equations (62) and (69) for two of the variables dt.» da, or 
db, setting the third equal to zero. It is possible to solve equations 


(62) and (69) for dt,, da, and db only if the rank of the matrix 


i lige 


0 
o4 a 
(70) "12 ) 
CBl "ama? $33 


is two, hence this is a criterion for deciding whether any admissible 
curves exist. If the rank of (70) is two, replace a, b, and a by a + da, 


b + db, and t, + dt, and repeat the calculations. Continue this process 


iL 
until either the rank of (70) becomes one or until some admissibility 
criterion is met. Note: In correcting two variables, it is essential 
that the value of third variable be chosen close to the value it will 
have on an extremal; otherwise the correction routine may not converge 
to an admissible path, even when such a path exists. The admissibility 
criterion used by the author was ESI - x(t,)| a \ * y(T)| § ens for 


the matrices used this criterion was normally met in less than ten iter- 


32 


ibs alia She 
ations when the determinants and went to zero, this 


%31 %35 %31 33 


was also obvious by the tenth iteration. 
When admissibility had been achieved, equations (62), (67), and 


(69) took on the form 


C71) 0 = a,,dt, + a19da 
(7 Deaax P= Qo dt, + Qooda a, .db 
(7 3) 0 = Qa, dt, + A25da + @32db. 


These equations were solved for dx(T) unless the rank of 


Q11 %12 Y 
| O11 %2 ° 
a Ao, 2 was equal to the rank of - If the rank of 
2 ee eS oe clmnne, 
Sle Cee 
O31 %32 %33 


the latter matrix was two, this condition was that the determinant of the 
first matrix must be zero. We therefore want some way to determine a new 
set of corrections to drive the determinant toward zero. The method 
adopted was to choose dt); da, and db proportional to the components of 
grad x(t,) X grad y(T), where "X" indicates the vector cross-proéuct and 
where the proportionality constant is some fixed number € multiplied by 
the determinant la; 5 | - The old values for ty.a, and b were replaced by 
the corrected values and a new admissible curve was found using the meth- 
od outlined above. For this curve the determinant lo; | and the correc- 
tions da, db, and dt, were again calculated together with x(T). This 
process of finding an admissible curve, then calculating a set of cor- 
rections to increase x(T), then finding a new admissible curve, then 
finding a new set of corrections to increase x(T) was continued until 

the determinant la; | changed sign. When the determinant changed sign, 


however, the above procedure no longer worked - the corrections became 


33 


larger rather than smaller, and the values of x(T) got smaller instead of 
larger. A procedure that worked when the determinant changed sign was 
the following: A new proportionality constant was chosen equal to one 
fifth of the sum of the positive and negative determinants; a new set of 
corrections da, db, and dt, was then calculated. The admissible curve 
calculated after finding this new set of corrections invariably yielded 
a larger x(T) than had been obtained before. This method of successive 
approximations seems to be a nearly foolproof method for determining the 
maximum x(T). The procedure was stopped when either the absolute value 
of the determinant was less than 107° or it became obvious that the con- 
vergence was too slow for the amount of time available on the computer. 
If convergence was too slow, it could usually be improved by increasing 
the number € referred to above. 

Another successive approximation procedure used was to compare the 
x(T) resulting from replacing t,,a, and b by their corrected values with 
the previous x(T). If the new x(T) were smaller, &was reduced by a fac- 
tor of ten and a new set of corrections and hence a new x(T) computed. 
This process was continued until either the sum of the absolute values 
of the corrections was less than 1075 or the new x(T) was larger than 
the old x(T), in which case the new x(T) became the value compared with. 
This process gives a monotonically increasing sequence of values of x(T); 
the speed of computation was comparable to that of the other method. 

It was necessary to insure that the corrections from the routine 
were not so large as to make the linearity approximations in the correc- 
tion integrals invalid. The author used the following procedure: In 
the correction to admissibility, if either da or db exceeded .3 in abso- 


lute value, or if dt, exceeded .5 in absolute value, he divided all the 
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~~ 


corrections in half. For the problems worked, this criterion was ade- 
quate although in several cases it slowed convergence quite a bit. (For 
example, where the initial value given for a was 1. and where the final 
value on the admissible curve was 12.8.) Once admissibility was achieved, 
it was necessary to use some linearity criterion on the second set of 
corrections. The procedure was to divide the corrections in half if the 
sum of the absolute values of the corrections exceeded .5. 

With the above correction procedures, computations were made on 
several sets of matrices with various boundary conditions. The only case 
where it was possible to verify the results by comparison with a physical 
situation where A and C were set equal to zero and where B and D were 


2 x 2 scalar matrices. This set of matrices gives the differential equa- 
‘ ie O<tct) 


tion X = an t]<teT’ The constraint on F was that Gp + C)- = 1; 


Pa and fo were consequently chosen as cos @ and sina, for 0* te t and 


12 
as cos £ and sin f®, for t,;<¢t<T. Solutions to the differential equation 
as set up describe the path taken by a light ray going from one isotropic 
medium into another, with @ and B being the angles between the light ray 
and the normals to the plane of discontinuity. As such, these solutions 
should obey Snell's Law, namely that the ratio of the velocities on oppo- 
site sides of the discontinuity plae should be the same as the ratio of 
the sines of the angles the ray makes with the normal to that plane. The 
cases tried were for Non = a x, = 1205 eos 5.0, a = 1.08andybp = 2.0 
and, in the second case, a = 2.0 and b= 1.0. For the first case, the 
ratio between the sines should be .5. & was found to be approximately 
-403378 radians and 8 approximately -902758 radians. Thecratio of the 


Sines of these angles was .500, which verified the accuracy of the rou- 


tine. For a = 2.0 and b = 1.0, the angles were reversed, and again the 
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accuracy of the routine was verified. | 


; a mss ap 2) 
Another case tried was where A= \, 4, 5}> B =| 9 1.3] » 


6 0 220 50.2 ; 
C= . on and D = 0 2.0}? with x, = 1.0, T = 4.0, and where imitéial 


values for t,, tana, and tan B were .73, 4.0, and .7, respectively. 

The routine converged to t, = -8935, tan @ = 1.4925, and tan 6 = .4967, 
with the maximum x(T) being 6.4622. This convergence took fifteen min- 
utes and fifty-one seconds on a CDC 1604 computer. The author felt that 
convergence could be improved if better initial estimates were made of 


t,, tana, and tan 8B. 
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APPENDIX | —, 
NOTES 
1.FLOW CHART PAGE NUMBERS ARE 


DESIGNATED BY THE LETTER F, eg., 
FG IS PAGE G OF THE FLOW CHART. 





2CONNECTORS ARE REFERENCED BY 
PAGE NUMSERS, @.9., re DESIGNATES 
CONNECTOR 3 IS ON PAGE F2. 






These constants are given: 
x(0), y(0), XMID, YFINAL, MATRICES A, 
B, C, and D, EPSFAC, X11, ATC12, ATD11, 
TONE, CGUESS, DGUESS, CAPT, ERROR, 
EPS1 







Calculate CKI = cy)” 


HANGES IN X(t,) AND Y DUE 
TO CHANGES IN CGUESS 

i CORR(1) 

XONE(2) = (1 0) CKI (ceee) 


YONE(2) = CORR(2) 


SET STEP SIZE 


TSTEP(1) = (CAPT - TONE)/50 
TSTEP(2) = TONE/50 







INTEGRATION FROM t = CAPT to t_= TONE 














VARIATION IN X(TONE) DUE TO A VARIATION 
IN TONE 


XONE(1) = (1 0) CKI 








PARTIAL VARIATION IN Y DUE 
TO A VARIATION IN TONE 


YONE(3) = CORR(4) 















6) 0, oF 


ee ae ee aw es ae ew ew es ow ww ew es awe ae es = 


Integrate: 












Calculate: 


SRT = ((acowi)? + (Acon@)) *) 


_ (ACON()) 
> Ctonea) 7 gs 






AX + BF 


Integrate: x 
d (corr(1)) = (-ccuzss x \Bi? x [«|? + (sRT)°) dt 






d (corR(2)) = (ja)? x |A|? + (SRT)*) at 






11 
21 


Kk! + ccuess x K’* 6x 


kK?! + ccuess x K72 6K 


-ccuEss x |p]? x |a| x DINTFAC + (SRT)? dt 


+ CGUESS x 6K?" 
+ CGUESS x 6K?" 








Calculate: DINTFAC = |B] 








Integrate: d (CORR(3)) 
d (CoRR(4)) = |B? x [A] x DINTFAC + (SRT)? at ke 





YDOT1 






x(t,) “ 
Out: ( ‘). K (t,.). Xit,) = (Ybor) CORRG), CORR), CORR(3), CORR(4) 


@n 


38 \F)) 


Yes 
[DELTAT| > 2? P| Set DELTAT = DELTAT/2 


NO 


# 


Print DELTAC, DELTAT 
Integrate: DELTAD , G} 


Calculate: (Acon(a) ACON(4)) 2 (xi? + DGUESS x «*7) D 


ip 


. »\1/2 
SRT = (aconis)) + (aconca)) 
TTWO = TONE + DELTAT 
_ faconis)\ . 
F = (ACON®?) ger 
Integrate: X = CX + DF = 

On first step only, set Se) = CX + DF TTWO <0? TONE = TONE/2 


Integrate: ; 2 2 3 ; 
D(CORA) =| -pGUESS x I|DI* x lxl° + (SRT)*) at 


D(CORB) =( [p|? x |Al? + (sRT)?}at 


; 


CGUESS = CGUESS + == 


DGUESS = DGUESS + -———, 





TTWO = CAPT ? me 


VARIATIONS IN X(CAPT) AND Y(CAPT) DUE 













TO A VARIATION IN TONE Yes 
1T,, , (XDOT1 - XDOT2 
XONE(5) = CORR(3) + K"7(¢,) (eer = Sacre ONE 
CGUESS = CGUESS + DELTAC 
XDOT1 - XDO DGUESS = DGUESS + 
YONE(1) = YONE(3) + K*7(t,) eben iets “ ae 
F3 
TONE =LONE+ CAPT (5 


2 







Changes in X(CAPT) and Y(CAPT) 
due to change in DGUESS are 


CGUESS = CGUESS + earn to x DELTAC 
CORB and CORA DELTAT 


DELTAT 
CAPT - TONE z DELTAD 


DGUESS = DGUESS + 2 DELTAT 


+ 





YDIFF 


XDIFF = XMID - XINMID | 
= YFINAL - YBND ©) a 






Yes 
~ [XDIFF] + |YDIFF] s ERROR ? _ Cs) J 
No 


it 


_ |XONE(1) YONE(1) 
DETER = |xoNnE(2) CoRRG) 


©: 


39 (F2) 


Print: TONE 
TOO LARGE 










Yes 
CAPT - TONE s .0001 ? 


No 


Calculate tan”! (CGUESS) and tan") (DGUESS) 


ATANC = ATANF (CGUESS) 
ATAND = ATANF (DGUESS) 


Print ATANC, ATAND 


Print: TONE 


Yes 
TONE < .0001 ? TOO SMALU 


No 
Print: ; 
CGUESS, TONE | 
lpevtcl + |pettp| <10°°? Yes. print: WE ARE DONE 
o: 
No 


(6 ) ATC1l. = ATANC 
: ATD11 = ATAND 


XONE (2): 0 XONE(]1) 













DELTC = ATANC - ATCll 
DELTD = ATAND - ATD11 





DTMT = |CORR(2) CORA YONE(1) 
CORR(1) CORB XONE(5) ey) = 





EPS = EPSFAC x DTMT 


¢ 


Yes 
x(T) > X11 ? x(T) is new X11 
Calculate DELT1T, DELT1C, DELTI1D So as to 
improve x(T) for admissible paths 


| 0 XONE(1) 
CORA YONE(1) 


4 XONE(1)  XONE(2) 
DELT1D = EPS x] YonE(1) — CORR(2) 


XONE(2) 0 
CORR(2) CORA | : 









¢ 





DELTI1C = EPS x 





CGUESS = CGUESS - DELTI1C 
DGUESS = DGUESS - DELTI1D 
TONE = TONE - DELTAT 





DELT1T = EPS x 








DELTIC] + [DELT1D| 
+ |DELTIT| < .5? 






No Divide all corrections 
by two ely F3 


Yes 


Print DELT1C, DELT1D, DELTI1T, EPS 


= CGUESS + DELTIC 
DGUESS + DELT1D & F3 


= FONE + DELTIT 





40 (F3) 





¥ 





Solve for DELTAC and DELTAD 


XDIFF\ _ /XONE(2) % DELTAC 
YDIFF CORR(2} CORA) \DELTAD 






No 
|DELTACI] < .3? 























Yes 
Solve for DELTAD and DELTAT aa ee 
XDIFF\ _ /XONE(1) 0 DELTAT 
YDIFF YONE()) CORA / \DELTAD 
©; 
y 
|DELTAD| < .3 ? IXONEW2) x CORRG)| < 10°59 Sarre) aa Yes CO) 
No 
y No 


Yes 


Divide all corrections by 
two 









Divide all corrections by 





Print: NO ADMISSIBLE 
PATHS EXIST 






two 





OL 


41 (F4) 


DAD 


1503 


501 


De i 


552 
502 
503 


iG: 


+AL= 
+, ,F5e2) 


APPENDIX II 


PROGRAM DISCON 
THIS PROGRAM COMPUTES MAX X(CAPT) FOR A GIVEN Y(CAPT) AND A GIVEN 
X(T1)sWHERE XDOT=AX+BF FOR T BETWEEN O AND T1sAND CX+DF FOR T 
BETWEEN Tl AND CAPT. CAPT 1S GIVEN,BUT T1 IS NOT. 
DIMENSION A292) 98292) 9C(292) »D( 292) sBK (292) sF (4) sTSTEP(2) sE (2925 
2)e9ELV(14) sDE(16) >» AAK(232)sDP(16) sACON( 4) sAKTB(292)5CORR 
(4) 9CKI(292)sFLES(2) »XONE(5) »sYONE( 4) »sBKA( 2592) »FMOR(2) sDELK(2 92) 


3 »DACON(4) »sDFE(4) »AK( 5516) 


READ 1503 ,EPSFAC 

FORMAT(1F20.10) 

READ 501 s( (A(T 9J) 9JF#192)s7T#=192) 9 ((B(I 9d) sJFl92) 5s T2192) 

FORMAT (8F10.1) 

READ 5015( (C(I 9J) 9J2192)9T=192) s((D( IJ) sJ=1l 92) sTHRl 92) 

PRINT 551 

FORMAT (14X91HAs29X91HB929X91HCs29X s1HD) 

PRINT 55250100 ACT 9J) s9J=1 92) 9 (BT 5J) sJE1 92) 9 (C(I 94) sJB1 92) (DI I 9d) sz 


+l52)sT=ls2) 


FORMAT (4(8XsF4e1s6X9F4e198X) ) 

READ 5025TONE »X0,YOsYFINAL »XMIDsCGUESS»,ERROR » 
FORMAT (8F10.8) 

READ 503,DGUESS 

FORMAT(F10.8) 

PRINT 5539TONE »X0sYOsYFINAL »XMID»,CGUESS »sERROR »CAPT 
FORMAT(4Xs12HTONE GUESS= 9F5e259 IQH(XOsYO)=(sF4elslHs sF4els 
9F4els7H XMID= 9F5e239H CGUESS= »F4e1s8H ERROR= 


CAPT 


SH) YFIN 
9—- 7¢696H CAPT= 


INTEGRATE ADJOINT BACKWARDS TO GET K(0) 
F(1l)=0. 

F(2)=0.5 

F(3)=0.5 

F(4)=1l. 


ATC11=56 
ATD11=56 


AND KINVERSE AT Tle 
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Xll=l.e 

EPS2=—-.01 
TSTEP(2)=TONE/50.6 
DO 1003 I=193 
YONE(1)=O6 
TSTEP(1)=(CAPT—-TONE) /506¢ 
TSTEP(1)=-TSTEP(1) 
TSTEP(2)=-TSTEP(2) 
BK(lsl)=le 
BK(291)=0.4 
BK(1s2)=064 

BK (2s2)=le 

T=CAPT 

DO 101 I=1l»52 

DO 101 J=l»s2 
E(2eTsJ)=A(I oJ) 
E(lselsJ)=C(IsJ) 
ELV(1)=BK (1591) 
ELV(2)=BK (1.592) 
ELV(3)=BK(251) 
ELV(4)=BK (292) 

DO 1004 IMA=558 
ELV(IMA)=0. 

IB=1 

DO 102 IE=1,5,50 

DO 103 IC=194 

DO 104 ID=154 


DE( ID) =ELV( ID) +F( IC) *AK(IC-1s5ID) 

DP (1)=-E(1Bs1s1)*DE(1)-E(1Bs2,1)*DE(3) 
DP (3)=-E(1Bs192)*DE(1)~-E(1Bs2,2)*DE( 3) 
OP (2)=-E(IBylsl)*DE(2)-E(1Bs2,1)*DE(4) 
DP(4)=-E(1B5y192)*DE(2)-E(I1Bs2,2)*DE(4) 


DO 103 ID=154 


AK(1CsID)=TSTEP(IB)*DP( ID) 


DO 105 1D=154 
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00 33 
0034 
0035 
0036 
O37 
00 38 
0039 
NN4NA 
NN41 
0942 
0043 
NN44 
0045 
0046 
ON47 
0048 
0049 
0050 
0051 
0052 
0053 
0054 
0055 
0056 
0057 
0058 
0059 
0060 
0061 
0N62 
0963 
0064 
N065 
0066 
0967 
0068 


105 ELV(ID)=ELV(C ID) +(AK(6 1910) +2¢*#AK(2510)+20¢*AK( 3510) +AK(451D) )/66 

102 T=T+TSTEP(IB) 
BK(1l»o1l)=ELV(1) 
BK (1s2)=ELV(2) 
BK (2s1)=ELV(3) 
BK (292) 8ELV(4) 
IB=2 

1002 ELV(5)=(A(191)—-Clls1))*¥BK(191)4+(Al291)-Cl(251))*BK(291) 
ELV(6)=(A(191)-Cl1l91))*BK(192)4+(A(291)-Cl291))*BK (292) 
ELV(7)=(A(6192)-C(192))*BK(191)4+(A(292)-C( 2592) )*BK (291) 
ELV(8)=(A(192)—-C(192))*BK(192)+(A(292)-C(292))*BK (292) 
DO 112 TA=15s50 
DO 903 IC=1>54 
DO 904 ID=158 

904 DEC ID)SELV(IO)+F( IC) *AK(IC-1510) 
DP(1)=-E(IBslsl)*OE(1)~-E(1Bs2,1)*DE(3) 
DP(3)=-E(I1Bs1ls2)*DE(1)~E(I1Bs2,2)*DE( 3) 
DP (2)=-E(IBsl191)*DE(2)—-E(1Bs2,1)*DE(4) 
DP (4)=-E(I1Bs192)*DE( 2) -E(1Bs2,2)*DE(4) 
DP(5)=-E(IBsls1)*DE(5)-E(1Bs2,1)*DE(7) 
OP (7)=-E(1Bs1s2)*DE(5)-E(1Bs2,2)*DE(7) 
DP (6)=-E(IBs191)*DE(6)—-E(1Bs2,1)*DE(8) 
OP (8)=-E(1B91s2)*DE(6)—-E(1Bs2,2) *DE(8) 
DO 903 ID=1,8 

903 AK(IC,ID)=TSTEP(IB)*DP(ID) 
DO 905 ID=198 

905 ELV(ID)=ELV( ID) +(AK( 1910) +2e*AK(291D)4+20*AK( 3,10) +AK(49I1D) )/66 

l12 T=T+TSTEP(IB) 

‘a BK(1ls1)=ELV(1) 
BK(1s2)=ELV(2) 
BK (2s1)=ELV(3) 
BK(2s2)=ELV(4) 
DELK(15s1)SELV(5) 
DELK(1s2)2ELV(6) 
DELK(2s1)2ELV(7) j 
DELK (252) =ELV(8) 
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DELK IS THE MATRIX OF DELK S AT T=0 
TSTEPCI=H—-TSTEP(1) 

TSTEP(2)=-TSTEP(2) 

THE K MATRIX» AS IS» IS AT T#FO, 

THE NEXT INTEGRATION IS FROM T=0O TO T=TONE 
ELV(1)=BK (191) 

ELV(2)=BK (152) 

ELV(3)=BK (251) 

ELV(4)=BK (292) 

ELV(5)=0.' 

ELV(6)=0. 

ELV(7)=X0 

ELV(8)=YO 

ELV(9)=DELK( 1591) 

ELV(10)=DELK(1s2) 

ELV(11)=DELK(291) 

ELV(12)=DELK (252) 

DO 202 IA=1;50 

DO 203 IC=1954 

DO 204 I[D=1914 

DE( ID) =ELV( ID) +F (IC) *AK(IC-15sID) 

DP (1)=-A(1591)*DE(1)-A(2591)*DE(3) 

DP (3)=—-A(1s2)*DE(1)-A(292)*DE(3) 

DP (2)=—-A(191)*DE(2)-A(2s1)*DE(4) 

DP (4)=-A(192)*DE(2)-A(292)*DE(4) 

ACON(1)= DE(1)*B(151)+DE(3)*B( 251) +CGUESS*(DE(2)*B(151)+DE(4)*B(2, 
*1)) 
ACON(2)=DE(1)*B(192)4+DE(3)*B(2,2)+CGUESS*(DE(2)*B(1s2)4+DE(4)*Bl 25 
+2)) 

SRT=SQRTF (ACON( 1) ¥ACON(1)+ACON(2)*ACON(2)) 
AKDET=DE(1)*DE(4)-DE(2)*DE(3) 
BDET=B(151)*B(292)-B(251)*B(1,2) 

DP (5)=-CGUESS*AKDET*AKDET*BDET#BDET/SRT¥*3 
DP (6)= AKDET*AKDET*BDET#BDET/SRT¥*3 
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Days 
9196 
ol Oy 
9108 
0199 
0119 


fo tail 


0112 
Ol 
0114 
Oii> 
0116 
0117 
0118 
O13 
0120 
Ola 
0122 
O123 
0124 
O22 
0126 
OL an, 
0128 
Olea, 
OL 2G 
gril ey: 
OV hey 
0133 
0134 
OJ 
0136 
OL3y7 
0138 
Ole 
0140 


— 


203 


206 


202 


DP (7)=A(191)*¥DE(7)4+A(192)*DE(8)+(B(151)*ACON(1)486192)*ACON(2))/SR 
+T 

DP (8)=A(251)*DE(7)4+A( 2592) ¥DE(8)+(8(251) *ACON(1)4+8(252)*ACON(2))/5SR 
+T 

DP (9)=-A(151)*DE(9)-A(291)*DE(11) 

DP (11)=-A(1s92)*DE(9)-A(2s2)*DE(11) 

DP (10)=-A(191)*DE(10)-A(291)*DE(12) 

DP (12)=-A(192)*DE(10)-A(292)*DE(12) 

DINTFAC=CGUESS* (DE (1) *DE(12)-DE (4) ¥DE(9)+DE(2)*0E(11)-DE(3)*DE(10) 


1 +CGUESS*(DE(2)*DE(12)-DE(4)*DE(10)) )-DE(1)*DE(11)+DE(3)*DE (9) 
DP (13) =-CGUESS*BDET*BDET#¥AKDET#¥DINTFAC/SRT¥*3 
DP (14)= BDET*BDET#AKDET*DINTFAC/ SRT ¥¥3 


DO 203 I1D=1514 

AK( IC, ID) =TSTEP(2)*DP( ID) 

DO 206 ID=1514 

ELV( ID) =ELVC ID) +°(AK(1sID)+20*AK( 2510) +20*AK( 3910) +AK(4910)) S66 
BK (191)=ELV(1) 

BK (ls2)=ELV(2) 

BK(291)=ELV(3) 

BK (292)=ELV(4) 

CORR(1)=ELV(5) 

CORR(2)=ELV(6) 

XMOD=ELV(7) 

YMOD=ELV(8) 

DELK(151)=ELV(9) 

DELK (152) =ELV(10) 

DELK(291)=ELV(11)— 

DELK (2592) =ELV(12) 

CORR(3)=ELV(13) 

CORR(4)=ELV(14) 

T=T+TSTEP(2) 
CORR 3 AND 4&4 GIVE THE VARIATIONS OF X AND Y AT T1 DUE TO DTi 
GET K INVERSE AT TONE AND F MINUS 
AKDET=BK(1591)*¥BK(292)-BK(1,2)*¥BK(251) 

XDOT1=DP(7) 

YDOT1=D0P(8) 
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559 





(304 


Motels 1 )=BKt252) /AKDET 

CKI(192)=-BK(291)/AKDET 

CKI(291)=-BK(192)/AKDET 

CKI(292)=BK(191)/AKDET 

FLES(1)=ACON(1)/SRT 

FLES(2)=ACON(2)/SRT 

YONE (3)=CORR (4) 
XONE(3)=CKI(191)*CORR(3)4+CKI(152)*CORR(4) 

XONE (1) =XONE(3)4+XDOT1 

XONE(2)= (CKI (151) *CORR(1)4+CKI(192)*CORR(2)) 
YONE(2)=CORR(2) 

XINMIOD=XMOD 

PRINT 5595XINMID 

FORMAT( 20H 559 
INTEGRATE FROM 
ELV(1)=BK (151) 
ELV(2)=BK (152) 
ELV(3)=BK( 291) 
ELV(4)=BK (292) 
ELV(5)=XMOD 
ELV(6)=YMOD 
ELV(7)=06 
ELV(8)=06 

DO 302 IA=1,50 
DO 303 IC=154 
DO 304 !D0=158 
DE(ID)=ELV( ID) 4+FC IC) ¥AK(IC-1ls1f) 

DP (1)=-C(191)*DE(1)-C( 291) *DE(3) 

DP (3)=—-C(1s2)*DE(1)-C( 292) *DE( 3) 

DP (2)=—-C(1ls1)*DE(2)—-C(291)*DE(4) 

OP (4)=-C(1s2)*DE(2)-C( 2592) *DE(4) 
ACON(3)=DE(1)*D0(151)4+DE(3)*0(2,1)+OGUESS*(DE(2)#0(151)+DE(4)*0( 251 


X AT TONE IS 
Tl TO CAPT 


»9F20-10) 


+)) 


ACON(4)=DE(1)*0(192)4+DE(3)*0(2,2)+DGUESS*(DE(2)*D0(1s2)+DE(4)*D( 252 


+)) , 


SRT=SQRTF (ACON( 3) *ACON (3) +ACON(4) ¥ACON(4) ) 
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all tee 
0178 
HLT9 
0180 
0181 
0182 
(o) ene: 
0184 
0185 
0186 
0187 
0188 
0189 
OLD 
0191 
O12 
0193 
0194 
0195 
0196 
O97 
Oh eye! 
Ole? 


" A2N0 


9201 
0202 
0203 
0204 
0205 
0206 
0207 
0208 
0209 
0210 
0211 
0212 


CKDET=DE(1)*DE(4)-DE(2)*DE(3) 
DDET=D(19s1)*D(252)-D( 192) *0(2,51) 
DP (5)=C(1l91)*¥DE(5)4+C(192) ¥DE(6)+(D( 191) *ACON(3)4D(192)*ACON(4))/5SR 
+7 
DP (6)=C(291)*DE(5)4+C( 292) *DE(6)+(D(291])*ACON(3)4+D(292) #ACON(4))/5SR 
+T 
DP(7)= DDET*DDET*CKDET*#CKDET/SRT ¥*3 
DP (8)=-DGUESS*DDET*DDET*CKDET#CKDET/SRT ¥*3 
DO 303 ID=158 
303 AK(IC,10)=TSTEP(1)*DP( ID) 
IF (IA-1)80158015802 
801 CONTINUE 
XDOT2=DP(5) 
YDOT2=DP(6) 
802 CONTINUE 
DO 306 ID=1,58 
306 ELV(ID) =ELV( ID) +(AK(191D) +20*AK(291D)+20*AK(391D) +AK(491D) )/60 
302 T=T+TSTEP(1) 
BKA(1s1)=ELV(1) 
BKA(152)=ELV(2) 
BKA(251)=ELV(3) 
BKA(292)=ELV(4) 
XMOD=ELV(5) 
YMOD=ELV(6) 
CORA=ELV(7) 
CORB=ELV(8) 
THIS COMPLETES THE INTEGRATION 
XONE(5)=CORR(3)+BK (191)*(XDOT1-XD0OT2)+BK(251)*(YDOT1-YDOT2) 
YONE(1)=CORR(4)+BK(152)*(XD0T1-]XDOT2)+BK (252) *(YDOT1-YDOT2) 
XBND=XMOD 
YBND=YMOD 
PRINT 5585XBND»YBND 
558 FORMAT(4H 558s9H X(T)= 9F20.1099H Y(T)® »4F20010) 
XDIFF=XMID-XINMID 
YDIFF=YFINAL=YBND 
IF (ABSF (XDIFF)+ABSF(YDIFF) -ERROR) 998399856001 
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5674 


9001 
3002 


' 
990 
562 


996 





DETER=XONE(1)*CORR(2)-XONE(2)*YONE(1) 
IF (ABSF (DETER) -1.2E-5) 600226002 56011 
DELTAC=XDIFF/XONE( 2) 
DELTAD=(YDIFF-CORR(2)*DELTAC)/CORA 
DELTAT=0-6 

IF (ABSF(DELTAC)~23)6671 96672 96672 
DEL TAC=e¢5*DELTAC 

DELTAD=.-5*DELTAD 

DELTAT=e¢5*DELTAT 

GO TO 6670 

GO TO 900] 

DELTAT=XDIFF/XONE(1) 
DELTAD=(YDIFF-YONE(1)*XDIFF/XONE(1))/CORA 
IF (ABSF (DELTAD)~-e3 )6674 3667536675 
DELTAD=e5*DELTAD 

DELTAC=-5*DELTAC 

DELTAT=-e5*DELTAT 

GO TO 6673 

CONTINUE 

DELTAC=0. 

DETER=XONE(1)*CORA ° 

IF (ABSF (DETER) -12E-5)919859198,9001 
IF (ABSF (DELTAT ) ~e5)990 9900239002 
DELTAT=e5*DELTAT . 
DELTAC=e5*DELTAC 

DELTAD=.5*DELTAD 

GO TO 9001 

PRINT 562sDELTACsDELTAT »sDELTAD 
FORMAT(12H DELTAC =) 9F20610914H DELTA Tl = 
= »9F 20210) 

TTWO=TONE+DELTAT 

IF (TTWO) 996 99969309 

TONE=TONE*e 5 
CGUESS=CGUESS+(TONE/DELTAT*DEL TAC) 
DGUESS=DGUESS+(TONE/DELTAT*DEL TAD) 
GO TO 310 


»F200105312H DELTAD 
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0249 
0250 
Campo) 
OZ 


BOZ2o¢ 


0254 
0255 
9256 
0257 
258 
9259 
0260 
OWS. 
0262 
0263 
0264 
0265 
0266 
N267 
268 
0269 
Q270 
A271 
O22 
O73 
0274 
O22 


0276 


0 2a 
0278 
Qiraa oo} 
0280 
0281 
0282 
0283 
0284 


309 IF(TTWO-CAPT 131129953995 
995 TONE=( TONE+CAPT)*-5 
CGUESS=CGUESS+( CAPT~TONE) *-5/DELTAT*DELTAC 
DGUESS=DGUE SS+( CAPT-TONE) *e5/DELTAT*OELTAD 
OG, 10 310 
311 TONE=TTWO 
CGUESS=CGUESS+DELTAC 
DGUESS=DGUESS+DELTAD 
310 CONTINUE 
IF (CAPT-TONE=~10E~-6)97715977199772 
9771 PRINT 9781 
9781 FORMAT(5H 9781%16H TONE TOO LARGE) 
GO T0"9991 
9772 IF (TONE~1.E-6)977399773 59774 
9773 PRINT 9782 
9782 FORMAT(5H 9782516H TONE TOO SMALL) 
GO TO 999] 
9774 CONTINUE 
322 CONTINUE 
PRINT 563sCGUESS»s TONE 
563 FORMAT (16H 563 NEW C IS) »F202¢10915H NEW TONE IS »5F20.10) 
999 CONTINUE = 
Go 190 1 
998 DTMT=XONE(2)* (CORA*XONE(5)-CORB*YONE(1) ) 
+ +XONE(1)*(CORB*¥CORR (2)-CORA*CORR (1) ) 
IF (ABSF (DTMT )—-201) 9601 5960139602 
9601 IF(ABSF(DTMT)—.00001) 960359603 59604 1 
9603 PRINT 9651 +XBND 
9651 FORMAT(5H 965199H XMAX= »9F20,10) 
GO TO 999] 
9604 CONTINUE 
GO TO 9602 
9602 EPS2=EPSFAC*DTMT 
DEL TAX=XBND-X11 
IF (DELTAX ) 8401284018402 
8401 EPS=EPS*.1 
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GO TO 9627 G32) 





402 EPS=EPS2 0322 
| X11=XBND 0323 
GO TO 9616 | 0324 

627 TONE=TONE-DELTI1T 0325 
: OGUESS=DGUESS-DELT1D 0326 
: CGUESS=CGUESS-DELT1C ps2 
9616 DELTIC=EPS*(=XONE(1)*CORA) 0328 
: DELT1D=EPS* (XONE(1)*CORR(2)—-XONE(2) *YONE(1)) 0329 
DELT1T=EPS* (XONE(2)*CORA) 0330 
5621 IF (ABSF (DELT1C)+ABSF (DELT1D) +ABSF (DELT1T)—05) 6622 56623 96623 0331 
$623 DELT1C=.5*DELTIC 0332 
| DELT1D=.5*DELT1D 0333 
: DELT1T=e5*DELTIT 0334 
: GO TO 6621 0335 
3? CONTINUE 0336 
PRINT 9611sDELT1CsDELT1DsDELT1T | 0337 

9611 FORMAT(5H 9601910H DELTAC= 95F202e10311H ODELTAD= ,4F20.10911H DEL 0338 
: MAT= @SFA0.20) 0339 
CGUESS=CGUESS+DELT1C . 0340 
DGUESS=DGUESS+DELT10 0341 
TONE=TONE+DELT1T 0342 
ATANC=ATANF (CGUESS) 0343 
ATAND=ATANF (DGUESS) 0344 
PRINT 9617sATANCsATAND 0345 

9617 FORMAT(5H 961736H C= 9F20¢10,10H RADIANSs s6H D= 45F20010s9H R 0346 
+ADIANS) 0347 
DELTC=ATANC-ATC11 0348 
DELTD=ATAND—ATD11 0349 

IF (ABSF (DELTC)+ABSF(DELTD)—-leE=6) 99813998159982 0350 

9981 PRINT 9983 . 0351 
9983 FORMAT(5H 9983949H CORRECTION IS LESS THAN 2000001 - WE ARE DON 0352 
1E//) . | 0353 

GO TO 9991 : 0354 
9982-2 ATC1L1=ATANC . 0355 
ATD11=ATAND 0356 


3) 


ae AS eS ~ a A A ES SS oa — A = A tert en ee, A RN 


9985 CONTINUE 
PRINT 9612sEPS 

9612 FORMAT(12H 9612 EPS=  »F20.10) 
AAK(191)=BK(191)+CGUESS*BK (192) 
AAK (291) =BK(291)+CGUESS*BK (292) 
AAK(192)=BK(191)+D0GUESS*BK (192) 
AAK (292) =BK(291)+DGUESS*BK (292) 
PRINT 9987 


9987 FORMAT(5H 9987918H K* AT TONE MINUS910Xe915HK* AT TONE PLUS//) 


PRINT 99889((AAK(I 9J) sJ=1l92) sTeloe2) 


9988 FORMAT(10X9F12e¢5914Xs9F12-5) 
PRINT 9989 


9989 FORMAT(1HO912H XDOT MINUS»10X,9HXDOT PLUS) 


PRINT 99909XDOT1sXDOT2sYDOT1sYDOT2 


9990 FORMAT (5X9F1l20e5910X9F1205) 
GO TO l 
9198 PRINT 9199 


9199 FORMAT(68H RANK OF FIRST CORRECTION MATRIX IS ONE - NO ADMISSIBLE 


»  LCURVES EXIST) 
9991 CONTINUE 
STOP 
END 
END 
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